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ABSTRACT 

We propose a method to estimate non-Gaussian error bars on the matter power spectrum 
from galaxy surveys in the presence of non-trivial survey selection functions. The estimators 
are often obtained from formalisms like FKP and PKL, which rely on the assumption that 
the underlying field is Gaussian. The Monte Carlo method is more accurate but involves the 
tedious process of running and cross-correlating a large number of N-body simulations, in 
which the survey volume is embedded. From 200 N-body simulations, we extract a non-linear 
covariance matrix as a function of two scales and of the angle between two Fourier modes. 
All the non-Gaussian features of that matrix are then simply parametrized in terms of a few 
fitting functions and Eigenvectors. We furthermore develop a fast and accurate strategy that 
combines our parameterization with a general galaxy survey selection function, and incor- 
porate non-Gaussian Poisson uncertainty. We describe how to incorporate these two distinct 
non-Gaussian contributions into a typical analysis pipeline, and apply our method with the se- 
lection function from the 2dFGRS. We find that the observed Fourier modes correlate at much 
larger scales than that predicted by both FKP formalism or by pure N-body simulation in a 
"top hat" selection function. In particular, the observed Fourier modes are already 50 per cent 
correlated at k ~ O.l&Mpc , and the non-Gaussian fractional variance on the power spectrum 
(cr 2 p /P 2 (k)) is about a factor of 3.0 larger than the FKP prescription. At k ~ 0.4/iMpc -1 , the 
deviations are an order of magnitude. 
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1 INTRODUCTION 

With new galaxy surveys probing a larger dynamical range of 
our Universe, our ability to constrain cosmological parameters 
is improving considerably. In particular, one of the most im- 
portant goal of modern cosmology i s to understand the nature 
of dark energy dAlbrecht et ail [2006), a challenging task since 
there are currently no avenues to observe it directly. It is how- 
ever possible to probe its dynamics via its equation of state to, 
which enters in the Friedmann equation that governs the expan- 
sion of the Universe. Among different ways a> can be measured, 
the detecti on of the baryonic acoustic oscillations (BAO) dila- 
tion scale (Eisenstein et al. 2005; Tegmark et al. 2006; Hiitsi 2006; 
IPercival et alj|2007l ; iBlake et alj|201ll) is one of the favorite, both 
because of the low systematic uncertainty and the potentially 
high statistics on can achieve with current {Huchra et alj 1 19901 ; 
lYork et alj|200ol ; IColless et al.ll2003l; iDrinkwater et alj|2010h and 
future galaxy surveys (Peterson et al. 2006; Acquaviva et al. 2008; 
Schlegel et all 120091; ILSST Science Collaborations et al.1 120091 ; 



Bemt ez et al. 2009 ; lBeaulieu et all20ld) . 
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The strength of the BAO technique relies on an accurate and 
precise measurement of the matter power spectrum, whose uncer- 
tainty i s propagated on to the dark energy parameters via a Fisher 
matrix (Teg ma^ll997h . It is thus of the utmost importance to have 
optimal estimators of both the mean and the uncertainty of the 
power spectrum to start with. The prescription to construct an es- 
timator for the power spectrum of a Gaussian random field, in a 
given galaxy survey, was p ioneered by Feldmann, Kaiser and Pea- 
cock (jpeldman e t all 19941) (FKP for short). It states that the survey 
selection function effectively couples Fourier bands that are other- 
wise independen t, and that t he underlying power should then be de- 
convolved dSato et alj201 lh . Th is technique has been used in power 
spectrum measuremen t such as {Feldman et alJl994l;[Percival et al.l 
1200 ll ; ICole et al.ll2005l ; lHutsill2006| ; |Blake et alj|20ld) . Althoughlt 
is fast, the error bars between the bands are correlated, plus it has 
the undesired tendency to smear out the underlying power spec- 
tram, which can effectively reduce the signal-to-noise ratio in a 
BAO measurement. In that sense, the FKP power spectrum is said 
to be suboptimal. 

The band correlation induced by the FKP prescription can 
be removed by an Eigenvector decomposition of the selec- 
tion function, fol lowing the Pseudo Karhunen-Loeve formalism 
dVogelev & Szalavlll996l) (PKL). This was used in the analysis of 
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the SDSS data dTegmark et alj|20o3) and is the most optimal (i.e. 
loss-less) estimator of a convolved Gaussian random field, as un- 
derstood from the information theory point of view. It is never- 
theless a well known fact that this Gaussian assumption about the 
field is only valid in the linear regime, since the non-linear gravi- 
tational collaps e of the density effectively couples different Fourie r 
modes together ( Meiksin & Whi te! 1999 l; lRimes & Hamiltonl20ol) 
and th e phases of the modes are no longer random l lColes & Chiang 
2000). Both the FKP and PKL prescriptions, by their Gaussian 
treatment, do not take into account the intrinsic non-linear coupling 
of the Fourier modes. It follows from this that for both methods, the 
measured power spectrum is suboptimal and the error bars are sys- 
tematically biased. Although the bias is usually small, it causes a 
problem when estimating derived quantities that need to be mea- 
sured with a percent level accuracy. 

For instance, the observed BAO signal sits right at the tran- 
sition between the linear and the non-linear regime, therefore the 
optimal estimator of the power spectrum must incorporate the 
non-linear modes. In particular, constraints on dark energy from 
BAO measurements require an accurate measurement of the mat- 
ter power spectrum covariance matrix. Under the FKP and PKL 
formalisms, the covariance matrix is biased as it tends to underes- 
timate the uncertainty and the amount of correlation between the 
power bands. Alternative ways of estimating the error, i.e. methods 
that involve mock catalogs, do model these non-linear dynamics, 
but it is not clear that the results are precise enough to measure 
four-points statistics, and we rather rely on accurate N-body simu- 
lations. 

Even more relevant is the recent realization that an optimal, 
i.e. non-Gaussian, estimate of the BAO dilation scale requires a pre- 
cise measurement of the inverse of the matrix, which is challeng- 
ing due the noisy nature of the forward matrix. It was nevertheless 
shown that, first, a suboptimal measurement of the power spectrum 
should be accompanied with error bars that treat the mean as subop- 
timal by consistency. These error bars diffe r from the naive Gaus- 
sian approximation by a significant amount iNga n et alj201lh . Sec- 
ond, it was shown in the same paper that an optimal measurement 
of the mean power spectrum could lead to an optimal measure- 
ment of the error, which differs from suboptimal measurements by 
a few percent. To achieve this, however, one needs to include the 
non-linear dynamics at a high precision, and needs in addition a 
strategy to cope with the noise present in the covariance matrix. 

There are thus two aspects of the analysis that need to be ex- 
tended from the above conclusions. On one hand, one would like to 
get rid of the aforementioned bias on the estimated uncertainty in 
actual data analyses, in which non-trivial survey selection function 
play a complex role. On the other hand, one would like to measure 
optimal error bars on the matter power spectrum of (non-Gaussian) 
galaxy fields. In this paper, we address the first issue, and provide 
a strategy to measure suboptimal but unbiased error bars on the 
power spectrum of galaxy surveys. It should be mentioned that our 
method could be put in conjunction with the PKL formalism to 
address the second issue, hence optimal measurement of both the 
power spectrum and of its error are now within reach. 

When constructing an estimator of the covariance matrix that 
corresponds to the sensitivity of a particular survey, the convolution 
with the survey selection function is one of the most challenging 
part. Whereas the convolution of the underlying power spectrum 
can be operated with angle averaged quantities, the convolution of 
the covariance matrix must be done in 6 dimensions, since the un- 
derlying covariance is not isotropic: Fourier modes with smaller 
angular separations are more correlated than those with larger an- 



gles dChiang et af]|2002l; iBernardeau et ai]|2002h . The first chal- 
lenge is to measure accurately this angular correlation, which is 
also scale dependent. Neither second order perturbation theory nor 
log-normal densities have been shown to calculate this quantity 
accurately, hence we must therefore rely on N-body simulations. 
This requires a special approach, since a naive pair counting of all 
Fourier modes in the four-point function, at a given angle, would 
take forever to compute. The second challenge comes from the 6- 
dimensional convolution of the covariance matrix with the survey 
function. This is a task that current computer clusters cannot solve 
by brute force, so we must find a way to use symmetries of the 
system and reduce the dimension of the integral. The fact is that 
the underlying covariance really depends only on three variables: 
two scales and the relative angles between the two Fourier modes. 
Moreover, it turns out, as we describe in section [6] that it is pos- 
sible to express this matrix into a set of multipoles, each of which 
can further be decomposed into a product of Eigenvectors. This ef- 
fectively factorizes the three dimensions of the covariance, hence 
the convolution can be broken down into smaller pieces. By doing 
so, the non-Gaussian calculation is within reach, and we present in 
this paper the first attempt at measuring deviations from Gaussian 
calculations, including both Poisson noise and a survey selection 
function. In short, the main ideas of this paper can be condensed as 
follow: 

(i) The underlying non-linear covariance matrix of the matter 
power spectrum exhibits many non-Gaussian features in the trans- 
and non-linear regimes. First, the diagonal elements of the angle- 
averaged covariance grow stronger, and correlation across different 
scales becomes important. Second, Fourier modes with similar (or 
identical) magnitudes correlate more if the angle between them is 
small. 

(ii) It is possible to model all of these non-Gaussian aspects 
(including the dependence on the angle between the two Fourier 
modes) with a small number of simple functions. 

(iii) With such a parameterization, it is possible, for the first 
time, to solve the six-dimensional integral that enters the convolu- 
tion of the covariance of the power spectrum with the galaxy survey 
selection function. 

Concerning the second point, the parameters that best fit our 
measurements are provided in section [7] but these are separately 
testable, and could be verified by other groups and with other ways. 
These are anyway expected to change when one uses haloes instead 
of particles. The third point is, however, a straight forward recipe 
that is robust under possible changes of best-fitting parameters, and 
provides, assuming that the input parameters are correct, an unbi- 
ased measurement of the non-Gaussian uncertainty of the matter 
power spectrum. 

As indicated by the title, this paper is the first part of a general 
strategy that aims at constructing an unbiased, non-Gaussian esti- 
mator of the uncertainty on the matter power spectrum measured 
in galaxy survey. The second part, which we thereafter refer to as 
HDP2 (in preparation), exploits the fact that the measurement of 
the C(k, k', 6) matrix provides a novel handle at measuring C(k, k')\ 
the two quantities are related by a straight forward integration over 
9. As shown in a later section of the current paper, it turns out that 
the main contribution to C(k,k') comes from small separation an- 
gles, while larger angles are noise dominated. It is thus possible to 
perform a noise weighted integral, which results in a more optimal 
measurement of C(k,k') and of its error bars, compared to direct 
or bootstrap sampling. It is therefore possible to extract an accurate 
non-Gaussian error bar on the power spectrum with a much smaller 
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number of realizations. This opens the door for a measurement of 
a non-Gaussian covariance matrix directly from the data (i.e. an in- 
ternal estimate), a significant step forward in the error analysis of 
galaxy surveys. In that second part, we will first test the improve- 
ment of the noise- weigthed technique on the same 200 simulations. 
We will then attempt to extract the mean and the error on all the el- 
ements of the C(k, k') matrix from a mock survey, that will consist 
of a handful of simulations. 

Back to the current paper, our first objective is thus to mea- 
sure the covariance of the power spectrum between various scales 
and angles, and organize this information into a compact matrix, 
C(k,k',6). We describe how we solve this problem in a fast way, 
which is based on a series of fast Fourier transforms and that can 
be run in parallel on a large number of computers. We found that 
the angular dependence, at fixed scales (k t kf), is rather smooth, 
it agrees with analytical predictions in the linear regime, but devi- 
ates importantly from Gaussianity for smaller scales. The depen- 
dence is somehow similar when the two scales are identical, up to 
an delta function for vanishing angles. We also found that, once 
projected on to a series of Legendre polynomials, it take very few 
multipoles to describe the complete original function. We perform 
this transform for all scale combinations and group the results in 
terms of the multipole moments. We also provide a recipe to recon- 
struct the original covariance matrix, with the angular dependence, 
from a handful of fitting functions. This is another advantage of the 
method we propose here: it can be separately tested and improved. 

Our second objective is to provide a general method to com- 
bine this C(k,k',0) with a survey selection function and non- 
Gaussian Poisson noise, and hence allow the extraction of non- 
Gaussian error bars on the measured power spectrum. We test 
our te chnique on the pub licly available 2dFGRS selection func- 
tions dNorberg et alj|20"02h and find that there is a significant de- 
parture between the Gaussian and non-Gaussian treatment. In par- 
ticular, the fractional error of the power spectrum (cr 2 p /P 2 (k)) at 
k ~ 0. l/iMpc~' is about a factor of 3.0 higher in the non-Gaussian 
analysis, and the departure reaches an order of magnitude by k ~ 
0.4/iMpc -1 . The method proposed here can be also applied to other 
kinds of BAO experiments, including intensity ma pping from the 
emission of the 2 1 cm line by neutral Hydrogen l lPeterson et al.l 
120061: lLazioll2008l: ISchlegel et alj|2009h. or Lvman-a forests sur- 
veys jMcDonald & Eisensteinl2007l;lMcOuinn & Whitel201 ll) . We 
did not, however, include the effect of redshift distortions, and fo- 
cused our efforts on dark matter density fields obtained from sim- 
ulated particles. An improved version of this work would include 
both of these effects, however. 

This paper is organized as follow: In section|2] we briefly re- 
view the FKP method, and describe how to estimate non-Gaussian 
error bars in realistic surveys, given a previous knowledge of 
C(k, k',9). We then lay down the mathematical formalism that de- 
scribes how we extract this quantity from simulated density fields 
in section [3] Section [4] contains the results from sanity checks and 
null tests that were performed to validate our method, and briefly 
describes our N-body simulations. We present our measurement of 
C(k, k', 6) in section[5] and describe the multipole decomposition in 
section|6] In section|7] we further simplify the results by extracting 
the principal Eigenvectors and provide fitting formulas to recon- 
struct easily the full covariance matrix. Section [8] contains results 
of applying our method for a set of simple selection functions. We 
finally discuss some implications and extensions of our methods in 
[9] and conclude in section [lOl 



2 MATTER POWER SPECTRUM FROM GALAXY 
SURVEYS 

In this section, we quickly review the g eneral FKP method , 
which is commonly used in data an alysis jFeldmanetai]|l994 
IPercival et alfeOoHlBlake et al.| [2010). We then point out some of 
the major the flaws of such techniques when measuring the un- 
certainty, and describe how non-Gaussian error bars could be es- 
timated in principle. Before moving on, though, we first lay down 
the conventions used throughout the paper. The reader familiar with 
the FKP method may skip to section [2~2l 

A continuous density field c5(x) is related to its Fourier trans- 
form c5(k) by 



S(k) 



(1) 



where k is the wave number corresponding to a given Fourier mode. 
The power spectrum P(k) of the field is denned as: 



<<5(k)<5*(k')> = (27r) 3 P(k)<5 D (k- k') 

and is related to the mass auto-correlation function by : 



(2tt) 3 



*) 3 J 



P(k)d 3 k 



(2) 



(3) 



In the above expressions, the angle brackets refer to a volume aver- 
age in Fourier space, and do(k) stands for the Dirac delta function. 

2.1 The optimal estimator of the power spectrum 

The power spectrum of the matter field contains a wealth of infor- 
mation about the cosmic history and the principal constituents of 
the Universe. Unfortunately, it is not directly detectable, since our 
observations are subject to cosmic variance, detection noise, light 
to mass bias, redshift distortions and incomplete sky surveys. The 
FKP method provides an optimal estimator of the matter power 
spectrum P(k) under the assumption that the density field is Gaus- 
sian. It is formulated in terms of the survey selection function IV(x), 
the galaxy number density n, the dimensions (n x , n y , w,) of the grid 
where the Fourier transforms are performed, and the actual number 
count per pixel n(x). All the following calculations can be found in 
dFeldman et alJI 19941) . and are included here for the sake of com- 
pleteness. 

The first step is to construct series of weights w(x) as 

vc (x) = = (4) 

1 + W(x)N c nP 1 + hP 

where N c - n x n y n z , h is the mean galaxy density and Pa is a char- 
acteristic amplitude of the power spectrum at the scale we want to 
measure. Since the latter is not known a priori, it is usually ob- 
tained from a theoretical model, and sometimes updated iteratively. 
The selection function is also normalized such that Yi\ W(x) = 1. 

The optimal estimator of the power spectrum, P„,(k), is ob- 
tained first by re-weighting each pixel by the weights in [Eq.[4), 
then by subtracting from the result a random catalog with the same 
selection function, weights and number of objects N. After taking 
the expectation value of the results, the 2-points statistics of the 
pixel counts becomes 



(n(x)n(x')) = hn'{\ + f(x - x')) + hS D (x - x') 



(5) 



where n is the mean density in the patch over which the average is 
performed. The Fourier transform is then given by 

|»(k) - NW(k)\ 2 - N W(x)w 2 (x) 



N 2 N c £ x W 2 (x)w 2 (x) 



(6) 
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where denominator is a convenient normalization. This measured 
power is aliased by the grid mass assi gnment sch eme, and should 
be divided by the appropriate function ( Jing 2005). 

What this estimator measures is not the underlying power 
spectrum P(k), but a convolution with the survey selection func- 
tion: 



<^.«Ck)> = 



I t , P(k')\W(k - k')| 2 
N c 2 X W 2 (x)w 2 (x) 



(7) 



It ideally needs to be deconvolved, an operation that is not always 
possible. 

For many survey geometries, the convolution effectively trans- 
fer power across diff erent bins which are uncoupled to start with 
(Tegm ark et al.l2006h . As mentioned previously, the PKL prescrip- 
tion also assumes that the density field is Gaussian, but rotates into 
a basis in which the bins are decoupled. In that sense, the PKL tech- 
nique is more optimal than the FKP, unless the selection function is 
close to a "top hat", in which case the induced mode coupling van- 
ishes. Both case, however, rely on the fundamental assumption that 
the underlying density field is Gaussian, which is known to be inac- 
curate in the trans- and non-linear regime, where one still wants an 
accurate measure of the power spectrum for a BAO analysis. Ob- 
taining accurate error bars is a requirement for optimal analyses, 
and we shall examine how these are usually obtained. 



2.2 The FKP covariance matrix 

The covariance matrix of the angle averaged power spectrum is a 
four point function that contains the information about the band 
error bars, and possible correlation between them. As mentioned 
earlier, it is required for many cosmological parameter studies. It is 
generally obtained from the power spectrum as 



C(k,k') = (AP(k)AP(k')) 



(8) 



where AP{k) refers to the fluctuations of the measured values about 
the mean, which is ideally obtained from averaging over many re- 
alizations. In a typical galaxy survey, such independent realizations 
are obtained by sampling well separated patches of the sky. Be- 
cause of the cost of such an operation, the number of patches is 
usually very small. The covariance matrix is thus not resolved from 
the data, and the error bars are obtained with external techniques, 
i.e. from mock catalog^ or directly from Gaussian statistics (see 
HDP2 for a prescription that overcomes this challenge). For a uni- 
form (top-hat) selection function, the Gaussian covariance matrix 
is estimated as: 



\k,k'). 



N(k) 



(P(k) + P shBt f6 b 



(9) 



where />,/,„, = 1 jn and N(k) is the number of Fourier modes that 
enters in the measurement of P(k). In the ideal scenario of perfect 

spherical symmetry and resolution, N(k) = 47rk 2 Ak[j^J , with Ak 
being the width of the fc-band. The Kronecker delta function en- 
sures that there is no correlation between different modes, an inher- 
ent property of Gaussian random fields. This equation can be easily 
be modified to deal with measurements without angle averaging. 

The FKP prescription provides a generalization of [Eq. |9) for 
the case where the selection function varies across the volume. It is 



1 We post-pone the discussion of mock catalogs until the next section 



obtained from [Eq.[6j and given by 
2 



C FKP (kX) 



N(k)N(k') 



-^IPGQf-lO + S'Ck-k')! 2 (10) 



where 



Q(k) 



£ x W 2 (x)w 2 (x)exp(/kx) 
£ x W 2 (x)w 2 (x) 



V(k) = , J_\Z* W(x)w 2 (x)exp(/kx) 



(11) 



(12) 

nN c j Zx W 2 (x)w 2 (x) 

In [Eq. [TO), P is taken to be the mean of the power spectrum at 
separation k - k'. Because the selection functions are usually quite 
compact about k = 0, that approximation is reasonable for Gaus- 
sian fields. Also, [Eq|9) can be recovered by setting W(x) = l/N c . 



2.3 Non-Gaussian error bars 

As mentioned in the last section, it is necessary to have access 
to many realization of the matter field in order to measure a non- 
Gaussian covariance matrix of power spectrum. This could in prin- 
ciple be done from data across many different patches in the sky, but 
even then, we have only one sky to resolve the largest modes, which 
would therefore be dominated by cosmic variance. Not to mention 
the cost and time involved in measuring many large but discon- 
nected volumes. Fortunately, N-body simulations are now accurate 
and fast enough to generate large numbers of measurements of the 
matter power spectrum. Since they model the non-linear dynam- 
ics of structure growth, the density fields they generate are non- 
Gaussian. The covariance matrix constructed from a high number 
of simulations indeed shows a correlation across different scales in 
the non-linear regime i Meiksin & Whi te 1999; Rimes & H amiltonl 
120051 iTakahashi et alj2009l : lNgan et alj|201 ll) . 

Although much more representative of the underlying covari- 
ance, such matrices are hard to incorporate in a data analysis, first 
because they are based on a fixed set of cosmological parameters, 
but also because the simulated volume is cubic and periodic. Each 
survey group typically needs to run at least one N-Body simulation, 
and measure the power spectrum with and without the measured se- 
lection function, in order to quantify the bias of their measurement. 
The complete approach would then be to run hundred of these to 
measure the covariance matrix, and to repeat for a range cosmolog- 
ical parameters values. This whole procedure is expensive, which 
explains why it is never done in practice. The alternative is to use 
mock galaxy catalogs, obtained, for example, from log normaliza- 
tion of Gaussian densities, second order perturbation theory (PT), 
haloPT, and so on. Unfortunately, the accuracy of such techniques 
at modeling the four-point functions and angle dependencies has 
not been fully quantified. 

Another artefact of the simulations is that the number of parti- 
cles can be arbitrary adjusted such as to suppress the Poisson noise 
down to a level where it is negligible. This is certainly not true for 
many galaxy survey, in which the number density is often much 
lower. We measure a non-Gaussian Poisson error by sampling ran- 
dom fields with a selection threshold chosen as to mimic the num- 
ber density of a realistic survey, and incorporate the effect manually 
in the analysis, as explained in section[8] 

To measure non-Gaussian error bars on a realistic survey, the 
most accurate procedure would be to convolve the best available 
estimator of the covariance matrix with the selection function. Be- 
cause the later is generally not spherically symmetric, it is the full 
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6-dimensional covariance matrix, C(k, k'), that needs to be inte- 
grated over. Let us suppose, for a moment, that we successfully 
measured that complete non-Gaussian covariance matrix. It would 
first contain an element for each Fourier modes k (i.e. with no an- 
gular averaging), and from [Eqs.[7]and[8), we can write: 



Cf\«(k, k ) : 



2 k „ k ,„<AP(k")AP(k'")W(k _ k")| 2 |W(k' - k'")P 



(N 2 N c ^W 2 (x)w 2 (x)) 2 



(13) 



where the angled bracket is nothing else but that full covariance 
matrix C(k",k"'). We can then simplify the result since the co- 
variance between two Fourier modes depends only the angle y be- 
tween them, but not on the absolute orientation of the pair in space. 
In other words, we make use of this symmetry argument to write 
C(k", k'") = C(k",k"', y) without lost of generality. This angle can 
further be expressed in terms of the two angles made by k" and k'" 
as 



cosy = cos0"cos0'" + sin$"sin$'"cos(0" - tp'") 



(14) 



We show in a later section of this paper that the true covariance 
matrix can be decomposed into a sum of factorized terms, each of 
the form F l {k")F 2 {k"')G x {ff' ,<p")G 2 {ff" So the double con- 
volution of [Eq. 1131 can actually be broken into a sum of smaller 
pieces, with at most 3-dimensional integrals to perform. 

This sets the path of this paper. We propose in the next section 
a novel technique to measure C(k" ,k"' , y), we next present the re- 
sults, and we finally perform the convolution on a realistic survey 
selection functions. 



3 MEASURING THE ANGULAR DEPENDENCE: THE 
METHOD 

As mentioned above, our first objective is to extract the covariance 
matrix of the power spectra from N-Body simulations, as a function 
of two scales and one angle: C(k, k',9). In this section, we develop 
a novel way to obtain covariances and cross-correlations and which 
allows us to perform this measurement. 

3.1 Cross-correlations from Fourier transforms 

We start by assuming we have measured the power spectrum from 
a large number of simulations. We first compute the mean of the 
angle averages: P(k) = (P(k)) NCi and the deviation from the mean 
of each mode: 



AP(k) = P(k) - P(k) 



(15) 



We then select two scales, k t and kj, that we want to cross-correlate. 
We make two identical copy of three-dimensional power spectra 
and multiply each one by a radial top hat function corresponding to 
the particular scales: 



AP,(k) = AP(kM|k|) 



(16) 



where Uj(k) = 6(k- ki)8(-k+ k t + 6k) is the product of two Heaviside 
functions. Also, 6k is the shell thickness, taken to be very small. We 
then cross-correlate the subsets and define: 



(17) 



E' J (Ak) = -i-r f AP,(k)APj(k + Ak)d 3 k 

We then express both AP,^(k)'s in lEq. 1171 in terms of their 
mass auto-correlation functions A£,j(x). We first integrate over 
exp[;k ■ (x + x')]d 3 k and obtain a delta function, which allows us to 



get rid of one of the real space integral. After slightly rearranging 
the terms, we obtain: 



£' 7 (Ak) = f A^(x)Af*(x)e"' Ak "'d 3 x 

In the above equation, A£ ; can be expressed as: 



(18) 



: AP(k)uj(\k\)d 3 k 



— tfdk e- ikx AP(k)dn 

n) 3 Jfc J 



(19) 



Since the shells we select are very thin, we can safely approximate 
that the power spectrum is constant over the infinitesimal range, 
and thus perform the k integral: 



A£(x) 



(2. 



X —£6k fe ik 

■n) 3 ' J 



s APi(k)dSl 



(20) 



We next repeat the same procedure for a scale j, multiply both 
auto-correlation functions together, and Fourier transform the prod- 
uct, following lEq. 1181 . The result is the cross-correlation S'-'(Ak), 
which becomes, after performing the x integral over the plane wave: 



(21) 
(22) 



E' 7 (Ak) = -— r*?ifc?<S 2 ifc fdnf dQ.'x 
(2ny J J 

APi(k)APj(k')6 D (k\ - ki - Ak) 

The delta function enforces Ak to point from kj to kj, which span 
the shells k t and kj respectively in the integrals over the solid angles. 
This geometry allows us to use the cosine law and relate |Ak| to the 
angle 6 it subtends, as seen in Fig.[T] such that: 



■ 



k) + k 2 - |Ak| 2 



(23) 



Since many Ak subtend the same angle 6, we can perform an aver- 
age over them and compute 

S''(0) = <S^(Ak))Ak=M (24) 



3.2 Normalization 

The quantity Y.''{9) is not exactly equal to C(k h kj, 9), because there 
is a subtle aliasing effect which is purely geometrical, and which 
needs to be canceled. To see how this arises, we work out a very 
simple scenario, in which the density field is perfectly isotropic. In 
that case, we can write AP(k) = AP(k), hence the angular integra- 
tion in rEq |20| is straight forward and we get: 

*? 

A?,(x) = A£(x) = -L-APMUkx) (25) 
txL 

with jo(x) the zeroth order spherical Bessel function. We have also 
assigned 6k = 2n/Lto the shell thickness, which corresponds to the 
resolution of a simulation of side L. Then, lEq |18l becomes 



where 

F ij (9) = J j Q (k,x)j {kjX)j {ex)x l dx 



(26) 



(27) 



The function F(k h kj,ff) is independent of the actual power spec- 
trum; it is purely a geometrical artifact, and we discuss in Appendix 
lAlits meaning. However, a power spectrum that is exactly isotropic 
should have no angular dependence.We thus define a normalization 
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Figure 1. Geometry of the system. For a fixed pair of shells, the magni- 
tudes of the Fourier modes k\ and ki are fixed, so the angle between them is 
directly found from the separation vector Ak. Note that we use interchange- 
ably number or roman letters to denote individual Fourier modes. 



3.3 Zero-lag point 

It is important to note that for a given realization, the point at = 0, 
which we refer to as the zero-lag point, must be treated with care. 
We recall from fEq. ! 191 that for a given simulation, we subtract the 
mean (P(k)) from every points which are on the specified /r-shell 
(and the other points are set to zero); we repeat the same procedure 
on a second fc-shell, and we then cross-correlate the two fields. 

When the two shells are identical, i.e. i = j, the zero-lag point 
of each simulation first computes the square of the deviation the 
mean P(k), then averages the result over the whole shell. It is equiv- 
alent to calculating the variance over the shell, but using a mean 
which is is somewhat off from the actual mean on that shell. That 
effectively boosts the variance. When we average over all simula- 
tions, the zero-lag points can be written as: 

E"(0) = (P 2 (k)) N , n - <P(k,)>^ Q (29) 

where, in the first term, the angle average and mean over all real- 
izations are computed after squaring each grid cell. By comparison, 
the variance on angle averaged power spectra would be obtained by 
performing, in the first term, the angle averaging first, then taking 
the square, then taking the mean. 

On the other hand, when the two shells are different, the zero- 
lag point is now the average over AP(k)AP(k') on both shells. Since 
we are no longer squaring each terms, it now includes negative 
terms, hence is generally of much smaller amplitude. 



E^(#), as the output of rEq. 1241 with AP(kjj) = 1 everywhere on 
the shells. The final results is obtained by dividing off this normal- 
ization, which cancels off the geometrical effect: 

C(k„ kj, 9) = ^ = (AP(k l )AP(k J )) (28) 

We stress again that this result is an average over all configurations 
satisfying kj = ki + Ak. To summarize, here is a condensed list of 
the steps taken to measure C(k, k', 6): 

(i) Measure the mean angle averaged P(k) from an ensemble of 
simulations, 

(ii) Select a combination of shells kij to cross-correlate, 

(iii) For each simulation, compute P(k), duplicate and multiply 
each replica by a top hat iiij(k), which effectively sets to zero every 
off-shell grid cells, 

(iv) Subtract P(k) from each cells in the shell, 

(v) Fourier transform both grids, complex multiply them, and 
Fourier transform back to fc-space, 

(vi) Loop over the Ak available, bin into E(|Ak| 2 ), and express 
the result as a function of 9, 

(vii) Repeat steps 5-6, but this time assigning the value of each 
cell in the shell to unity. Divide 1.(9) by this normalization. This is 
a measure of C(fc, , kj, 9) from one simulation, 

(viii) Repeat for all simulations, then compute the mean, 

(ix) Iterate over steps 2-8 for other shell combinations. 

To achieve better results, we make use of the fact that P(-k) = 
P(k), hence, following [Eq[T7), we can write E ij '(-Ak) = E'-'(Ak). 
This translates into a theoretical symmetry about 9 = n/2 in the 
angular dependence of the covariance. That property turns out to be 
very useful at reducing the numerical noise, since we can measure 
the covariance over the full angular range, but fold the results on to 
< 9 < k/2. Also, to avoid interpolating error, we chose to bin in 
(Ak) 2 before transforming to 9. 



4 VALIDATION OF THE METHOD 

We describe in this section a series of tests that compare our numer- 
ical results to semi-analytical solutions. We apply our numerical 
methods on a few simple situations in which we control either the 
density field or the three dimensional power spectrum. We first test 
our recipe on a power spectrum that is set to the second Legendre 
polynomial. The outcome can be compared to semi-analytical cal- 
culations and gives a good grip on the precision we can achieve. We 
next measure the angular dependence of the covariance matrix of 
white noise densities by Poisson sampling many random distribu- 
tions, and_present an estimate for a non-Gaussian Poisson error See 
dCohrJ 2006) for discussions of these two types of noise in a cos- 
mological context. We finally measure the angular cross-correlation 
from Gaussian random fields in order to later measure departures 
from Gaussianity. 



4.1 Testing C(k, k',9) with a Legendre polynomial 

As a first test, we enforce the z-dependence of the power spectrum 
to be equal to the second Legendre polynomial, and then compare 
our results to semi-analytic predictions. We manually set P(k) = 
kl, which is thus constant across the x — y plane. The mean and 
the deviation from the mean on a shell are given by {P (k))n = 
k 2 /3, AP(k) = (2/3)k 2 P 2 Qi) respectively, where P t (x) is the f-th 
Legendre polynomial and /j is the cosine of the inclination angle. 
The mass auto-correlation function associated with this power is 

-2k 4 

= -rrh(k,x) (30) 

The angular dependence of the covariance can be calculated semi- 
analytically from Eq[T8]and the above mass auto-correlation func- 
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Figure 2. (top:) Angular dependence of the covariance of a power spectrum 
set to the 2 Legendre polynomial, calculated here at k l= j = l.OAMpc . 
The solid line is the semi-analytical prediction. The curve is normalized 
to the value of the zero-lag point, such that it represents the actual cross- 
correlation coefficient between the Fourier modes. In this case, modes that 
points in like-directions are strongly correlated, (middle :) Angular depen- 
dence of the power spectrum cross-correlation coefficient measured from 
200 Poisson sampled random fields. The error bars were obtained by 500 
bootstrap sampling. We have selected two ^-shells i, j that are off by one 
grid cell: kj = k, ■ + 6k, with 8k = 0.0314/i/Mpc and it, ~ 1.0/iMpc -1 . 
The distribution for (i = j) is similar in shape, except for the zero-lag 
point which is much larger than any other points, and the plateau that is 
slightly higher. The solid line in this figure is the predicted value, which 
is well within the error bars. We have reproduced a similar plot for Pois- 
son densities with 8.0 million peaks, which is also flat, and find that the 
height of the plateau scales roughly as 1 /n 3 , where n is the number of Pois- 
son sampled objects, (bottom:) Angular dependence of the power spectrum 
cross-correlation coefficient, measured from 200 Gaussian random fields, 
this time with k: = k, + 5Sk , and again k{ ~ l.OAMpc . The theoretical 
prediction is zero, whereas we measure a constant 6 per cent correlation 
bias across all angles. We have verified that this bias is scale independent 
by changing ku. 



tion. The angular integration is straight forward, and we obtain 



Z, ij (Ak) : 



9nL Jo 



(3D 



We perform the x integral with k i= j = l.O/zMpc" , repeat 
the procedure for X^(Afc), and obtain a semi-analytical prediction: 
C(k, k' ,8) ~ _Pi(cos0), up to numerical noise. This agrees well with 
the numerical results produced by our technique, as shown in the 
top part of Fig. [2] We are plotting the angle-dependence of the co- 
variance matrix, normalized by the angle average of the covariance, 
such that the curve represents the actual cross-correlation coeffi- 
cient between the Fourier modes. We mention here that in the case 
where kj t kj, which we encounter in the following sections, we 
normalize to the square root of the product of the corresponding 
matrix elements: 



r(ki,kj,0). 



C(k„kj,8) 



T/C(k,,ki)C(kj,kj) 



(32) 



4.2 Testing C(k, k',9) with Poisson-sampled random fields 

To measure the response of our code to white noise, we produce a 
set of 200 density field representing Poisson sampling of random 
distributions. The sensitivity threshold is chosen such that ~ 8000 
peaks were counted on average. The standard deviation in the mea- 
sured P(k) decreases roughly as k~ 2 , expected from the fact that the 
number of cells on a fc-shell grows as k 2 . 

Because of the random nature of Poisson density, the variance 
on a given shell should be roughly constant across all directions. 
Moreover, after averaging over many realization, Poisson densi- 
ties are in principle statistically isotropic. We thus expect the mea- 
sured angular dependence of the covariance to be very close to flat, 
and, from [Eq |28l . we estimate it should plateau at a value some- 
what similar to the covariance of angle averaged power spectrum, 
C(k,k'): 



Cpoi SS onik,k',fj) ~ Cp oisson (k,k') + AStk'S, 



y,±i 



(33) 



In the particular case under study in this section, the Fourier modes 
separated by small angles are strongly correlated by construction. 



where yu = cos# and the two delta functions ensure that modes 
with different directions or scales do not couple together. Also, 
Cpoisson(k, kf) is obtained directly from the covariance matrix of the 
angle average power spectra. The constant A is much larger than 
Cp i S son(k,k'), for reasons explained in section |33l but the precise 
value is irrelevant to the current analysis. Fig. [3] shows the cross- 
correlation coefficient matrix for non-Gaussian Poisson noise. We 
observed that the angle-averaged modes are correlated by more 
than 30 per cent between scales smaller than k = 1.0/i/box. The 
reason for this feature is actually independent of cosmology, even 
though the matrix has a look very similar to that measured from 
simulation^] . The explanation lies in the fact that our Poisson den- 
sities do not have exactly the same number of objects, hence the 
asymptotic value of P(k) is not a perfect match for all field. This 
slight scatter in power translates into a correlation between the high 
fc-modes of a given density field, and we use this matrix to estimate 
the non-Gauss ian Poisson n oise. This is in good agreement with the 
predictions of dCohnl2006h . which calculated that the Poisson sam- 
pling of Gaussian fields induce non-Gaussian statistics, and that 
well separated scales can correlate significantly. 

We then measure the angular dependence of the covariance 
for these 200 Poisson distributions, also at k ~ 1.0/zMpc -1 . We ob- 
tain a distribution which is indeed close to flat, and consistent with 
a uniform 10 per cent correlation, as shown in the middle plot of 
Fig. [2] As before, we have normalized the plot by the angle av- 
eraged covariance matrix element, such as to exhibit the angular 
cross-correlation. Because the zero-lag point is typically a few or- 
ders of magnitude above the other points, we quote its value in the 
text or in the figures' caption where relevant, and resolve the struc- 
ture of the other angles. The mean of the un-normalized distribution 
is 133.3Mpc 6 /r 6 , a 10 per cent agreement with our rough estimate. 
We have re-binned the distributions on to a set of points that are op- 
timal for the upcoming angular integration, as described in section 
i 



4.3 Testing C(k, k', 6) with Gaussian random fields 

The next test consists in measuring the angular dependence of the 
covariance from of 200 Gaussian random fields. We use 200 power 



2 It is in fact arguable that such a matrix, constructed from a set of Pois- 
son densities, could have better performances at modeling the 'true' non- 
Gaussian covariance matrix, compared to the naive Gaussian approxima- 
tion. 
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Figure 3. Cross-correlation coefficient matrix, measured from the power 
spectra of 200 Poisson sampled random densities fields, selected to have 
8000 peaks on average. The correlation in high &-modes is unphysical, as 
explained in the text, and this represents our estimate of the non-Gaussian 
Poisson uncertainty. 
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Figure 4. (top:) Power spectrum of 200 simulations, produced by cubep3m, 
compared to CAMB at z = 0.5 (solid line). The error bars are the l<x stan- 
dard deviation on the 200 measured P(k). We only include modes with 
k < 2.34fcMpc in this analysis, as indicated by the arrow in the figure. 
(bottom:) Ratio between the simulated and modeled power spectra. 



spectra measured at z = 0.5, obtained from N-Body simulations 
(section [5~TT ), to generate 200 fields. Similarly to the Poisson fields, 
we expect the distribution to be be overall flat, except for the zero- 
lag point. Because we choose not to Poisson sample these Gaus- 
sian densities, the randomness should be such that near to perfect 
cancellation occurs between the different angles, and the plateau 
should be at zero. In the continuous case, the Gaussian covariance 
can be expressed as [Eq. 1411 . with 

C G m, ss (k h kj,n) = ^ ' SijSp^i (34) 
N(ki) 

where N(k) is the number of Fourier modes in the fc-shell. For k t = 
kj, the zero-lag point contains perfectly correlated power, so we 
expect it to have a very large value. As explained in section [3~3l we 
cannot directly compare its value to 2P 2 (k)/N(k), since the former 
is bin dependent, while the latter is not. In the case where i + j 
however, the zero-lag point should drop down to numerical noise. 

The measured angular dependence is presented in the bottom 
part of Fig. [2] where we see that the distribution is flat and con- 
sistent with 6 per cent correlation. This indicates that our method 
suffers from a small systematic bias and detects a small amount of 
correlation, in a angle independent manner. We have repeated this 
measurement for different scales ktj and obtained the same bias. 
We therefore conclude that any signal which is smaller than this 
amount is bias dominated and not well resolved. 



5.1 N-body simulations 

Since our Universe is not Gaussian at all scales relevant for BAO or 
weak lensing analyses, a robust error analysis should be based on 
non-Gaussian statistics, and, as mentioned earlier, N-body simula- 
tions are well suited to measure covariance matrices. Our numerical 
method is fast enough that, for fixed kj and kj, we can compute the 
angular dependence of the covariance matrix in about one minute. 
The average over 200 realizations can be done in parallel, hence 

producing all available combinations takes very little time. 

The simulations are produced by cubep3m( lMerz et aill2005h . 
a public N-body code that is both openmp and mpi parallel, which 
makes it among the fastest on the markeJJ We generate 200 Gaus- 
sian distributions of 200 Mpc/T 1 per side, with 256 3 particles, 
starting at z, = 40, and evolved them until z = 0.5. The simu- 
lations are run on the CITA Sunnyvale cluster, a Beowulf clus- 
ter of 200 Dell PE1950 compute nodes, each equipped with 2 
quad cores Intel(R) Xeon(R) E5310 @ 1.60GHz processors. Each 
node has access to 4GB of RAM and 2 gigE network interfaces. 
The power spectrum of these simulations is shown in Fig. [4] and 
shows a good agreement wit h the non-linear predictions from CAMB 
dSeliak & Zaldarriaga|[l99^) , up to k ~ 0.25/iMpc -1 . Beyond that 
scales the structures are under estimated due to the resolution 
limit of the simulations. For the rest of this paper, we only con- 
sider well resolved scales, in occurrence those in the range k e 
[0.314, 2.34]frMpc~', which we organize into 75 linearly spaced 
bins. 



5 MEASURING THE ANGULAR DEPENDENCE 

In this section, we present the results obtained from the measure- 
ment of the angular covariance in our 200 simulations. We explore 
different scale combinations and attempt to compare the outcome to 
expected results whenever possible. In particular, the linear regime 
should somewhat reproduce the behavior we observe in Gaussian 
random fields (see section |4~3l . In all figures, the error bars were 
obtained from 500 bootstrap re-sampling of our simulations, un- 
less otherwise specified. 



5.2 Results 

We present in Figs.[5]and[6]the angular dependence of the covari- 
ance between the power spectrum of various scales. As explained 
in the previous section, the distributions are normalized such as to 
represent the cross-correlation coefficient between modes separated 



3 http://www.cita.utoronto.ca/mediawiki/index.php/CubePM 
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Figure 5. Angular dependence of the power spectrum cross-correlation, 
measured from of 200 density fields, at &, =; = 0.17,0.46,0.93 and 
2.34/iMpc~' . The distribution exhibits a correlation of less than 10 per cent 
for angles larger than about 30". For scales less than 0.5/iMpc~' , the corre- 
lation increases up to 15 per cent for angles smaller than 10°, and to more 
than 40 per cent for smaller angle. 



Figure 6. Angular dependence of the power spectrum cross-correlation, 
measured from of 200 density fields, at k t =0.61, and kj = 0. 14, 0.46, 0.93 
and 2.34/iMpc~ 1 . The distribution exhibits a correlation of less than 10 per 
cent for angles larger than about 30° . For scales of similar sizes, the corre- 
lation increases up to 15 - 20 per cent for angles smaller than 15°. 



by an angle 6. In the first figure, both scales are selected to be iden- 
tical, and vary progressively from k = 0.17/iMpc~' to 2.34/iMpc~'. 
Modes separated by an angle larger than 30° are less correlated at 
all scales, and the correlation is even smaller for modes smaller 
than 0.5/iMpc~'. These latter modes are grouped in larger bins due 
to the higher discretization of the shells, and ideally one would like 
to run another set of simulation at larger scales to have a better reso- 
lution on those scales. However, the modeling of these rather larger 
scales have very little impact on the non-Gaussian analysis we are 
carrying, we therefore do not attempt to improve the situation. For 
highly non-linear scales, the correlation between modes less than 
10° increases up to 55 per cent. 

In the second figure, one of the two scale is held constant, 
at k = 0.61/iMpc~', while the other varies over the same range. 
Modes separated by angles larger than 30° are less than 10 per cent 
correlated, for all combinations of scales. When the two scales are 
of comparable size, the the correlation climbs up to values between 
15 and 20 per cent for angles smaller than 15°. 

This angular behavior is enlightening, as it shows how the er- 
ror between Fourier modes separated only by a small angle tend 
to correlate first. Qualitatively, this validate the fact that in non- 
Gaussian densities, quasi-parallel Fourier modes are probing essen- 
tially the same collapsed structures. When the angle is closer to 90°, 
however, one mode could go along a filament and the other across 
it, producing only weak correlations. It could thus be possible to 
construct a highly clustered density in which we could observe an 
anti-correlation at 90°, provided we are not noise dominated. 

This coherent behavior is a clear sign that the non-linear 
structure underwent gravitational collapse, and the departure from 
Gaussianity and white noise is obvious. Another signature of non- 
Gaussianity is that even in the presence of a small offset between 
the scales, the small angle correlation has a value higher than those 
at larger angles, because of the coupling between those scales. Fig. 
[6]shows this effect. 



5.3 From C(k h kj, 6) to C(k h kj) 

It is possible to recover the covariance matrix one obtains from 
the angle averaged P{k) by performing a weighted sum over the 
angular covariancfl Another test of the accuracy of our method is 
thus to compare the two ways that measure C(k,, kj). This is by far 
the least convenient way of measuring this matrix, and we perform 
this check solely for verification purposes. 

We perform this weighted sum and construct C(ki,kj), then 
compute a similar matrix from our 200 angle averaged power spec- 
tra. We present in Fig.|7]the cross-correlation coefficient matrix (see 
lEq. 1321 ) obtained in the first way, and show the fractional error be- 
tween both methods in Fig. [8] We observe that they agree at the few 
percent level, so long as we are in the non-linear regime. At very 
low fc-modes, however, many matrix elements are noisy due to the 
discretization of the shell; the (Ak,8) mapping in this very coarse 
grid environment becomes unreliable, and the re-weighting hard to 
do correctly. This results in high fractional errors, but at the same 
time, this region is still in the regime where the analytic Gaussian 
prediction is valid. In addition, this paper attempts to solve the bias 
caused by the non-Gaussianities that lie in the trans-linear and non- 
linear regime, in which discretization effects are much smaller. Fi- 
nally, we recall that these matrix elements have very little impact on 
most parame ter studies since such scales contain almos t no Fisher 
information teimes & Ha milton 2005l lNgan e"tal]|201lh . 

6 MULTIPOLE DECOMPOSITION 

As shown in last section, we have extracted the covariance ma- 
trix C(kj,kj,8) of matter power spectrum, cross-correlating the 75 

4 The weight here is simply the number of contribution that enter each an- 
gular bin, divided by the square of the total number of cells on the fc-shell. 
In other words, because the angular covariance we measure is an average 
over many pairs of cells, that average must first be undone, then the dif- 
ferent angles are summed up, and we finally divide by the total number of 
contributions. 
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Figure 7. Cross-correlation coefficient matrix, as measured from the angu- 
lar covariance. Each matrix element i, j was obtained from a reweigthed 
sum over C(kj,kj,6). T his is consistent with m a trices previously mea- 
sured in the literature fames & Hamiltonl |2005|; iTakahashi etafl |2009| ; 
iNgan et aljioIH) 
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Figure 8. Fractional error between the covariance matrices obtained with 
the two methods. We have suppressed the largest scales, which are noisy due 
to low statistics, and present the percent level agreement at smaller scales. 



different scales selected. Since the final objective is to incorporate 
this massive object into generic data analysis pipelines, it must be 
somehow simplified or made more compact. A quick glance at the 
figures of section [5] reveals that the angular dependence of the co- 
variance can be decomposed into a series of Legendre polynomials, 
in which only a few multipoles will bear a significant contribution. 
This allows us to rank the multipoles by importance and to keep 
only the dominant ones. These results are further simplified in sec- 
tion|7] where we provide fitting formulas to reconstruct C(kj, kj, 9). 

In this section, we describe how we perform this spherical har- 
monic decomposition, then we test our method on the control sam- 
ples described in section [4] and we finally measure the C({k,k') 
from the simulations. 



6.1 From C(k h kj, 9) to C t (k h kj) 

Here we lay down the mathematical relation between C(kj,kj,9) 
and C((kj,kj). Let us first recall that the spherical harmonics 
Y {m (G, (/>) can serve to project any function F(9, if>) on to a set of 
a Cm as: 



atm = J Y ( '" 



(9, <p)F(9, (/>)dCl 



(35) 



We substitute F(£2) -> AP,(k) = AP t (k, fl), which causes the coef- 
ficients to be scale dependent, i.e. at m — > a{ m (k). The angular power 
spectrum at a given angular size 9 ~ l/£ is defined as 



1 i 

C ( (kj,kj) = - 2_j \ a em{ki)a* e Jkj)\ 



(36) 



Combining both equations, and writing C'J = C((Jq,kj) to clarify 
the notation, we get 



q = — !— 



^ j Y""'(£l')Y e '"(Q.)x 



m=—€ * 

AP(k h £l)AP'(kj, Q!)dQ.(Ki 



(37) 



We use the completion rule on spherical harmonics to perform the 
sum: 



^ Y em (£l)Y <m (£l') = Hlip^cosy, 



An 



(38) 



(40) 



(41) 



where y is the angle between the £1 and fl' directions, and where 
P((x) are the Legendre polynomials of degree I. We then write 

C'J = J AP(k i ,£i)AP , (kj,n')P c (cosy)d£lclQ.' (39) 

Since we know that ki + Ak = kj, we make a change of variable 
and rotate the prime coordinate system such that k always points 
towards the z-axis. In this new frame, we have dQ." = dcos9"dif>", 
where 9" is the angle subtended by Ak. 9" thus corresponds to the 
angle between the two Fourier modes k and k'. It is also equal to y 
in rEq. |38| .We perform the 'unprime' integral first, which gives 

C'/ = -L T P t (cosy) J APj(k)APj(k + Ak)dndQ." 

The inner integral is C(k h kj, y), we rename y — > 6 and obtain 

Cf = J P e (casQ)C{k h kj,ff}dQ. 

In practice we are dealing with a discretized grid, hence we 
must convert the integral of lEq |41l into a sum. To minimize the 
error, we use a Legendre-Gauss weighted sum, the details of which 
can be found in the Appendix. In order to validate our method, we 
designed a few tests which are explained in the following sections. 



6.2 Testing Cc with a Legendre polynomial, with Poisson and 
Gaussian distributions 

We start our tests by measuring the Cf(k h kj) from the angular de- 
pendence of the covariance of power spectra, which is explicitely 
set to the second Legendre polynomial on the selected ^-shells, as 
described in section |4~T1 We expect the projection to produce a delta 
function at I = 2, up to numerical precision, since the Legendre 
polynomials are mutually orthogonal. We observe from this simple 
test a sharp peak at t = 2, which is about two orders of magnitude 
higher than any other points. 
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a slight offset between the two scales, in which case the analytical 
prediction is zero. Our results are consistent with zero for all multi- 
poles except I = 0, which receives an undesired contribution from 
the constant 6 per cent bias described in section l4~3l and observed in 
Fig-IH It turns out that this Co contribution is very small compared 
(i.e. less than one per cent) to the values obtained from simulated 
density fields, hence we do not attempt to correct for it. In the case 
where the two shells are identical, we observe similar results, up to 
an upward shift caused by the zero-lag point, which propagates to 
all multipoles. 

When performing this decomposition on the covariance ma- 
trix obtained from actual density fields, the departure from Gaus- 
sianity can be quantified as the number of distinct t needed to de- 
scribe the angular power distribution. 



Figure 9. (top :) Angular power of the cross-correlation obtained from 
200 Poisson densities, at kj^j ~ l.OAMpc , with an offset of one grid cell 
between the two scales, corresponding to 6k = 0.0314/iMpc~'. The power 
at t + is consistent with zero, as expected from [Eq. 1421 . We recall that 
the angular dependence of the covariance from Poisson densities is very 
weak, hence it projects almost exclusively on the ( = term, (bottom :) 
Gaussian angular power at k\ ~ 1.0/iMpc -1 , and^ = k, +56k. The analytical 
prediction is zero at all multipole, while we measure a Co term of about 
80.5/i~ 6 Mpc 6 . This is caused by the 6 percent bias we observed in Fig.ff] 



We next measure the Cc from the covariance matrix of Poisson 
densities, whose angular dependence, we recall, is close to flat (see 
section l4~2t , up to a delta function on the zero-lag point when the 
two shells are identical. From the orthogonality of the Legendre 
polynomials, a flat distribution is projected exlusively on the first 
multipole, we thus expect C Pmsso "(k t k') to peak at I = 0, and 
to vanish for other I. Moreover, we expect the C Pmsson (k = k') to 
exhibit, in addition, a vertical shift caused by the integration over 
the zero-lag point. The analytical expression can be obtained from 
rEqs. I33I4T1 . The azimutal integration gives a factor of 2n, the \i 
delta function gets rid of the last integral, and we get: 



C Pmssm (k,k') = 2nC Pm 



,ki=k' 



C e °' ss °"(k,k') = InCpaus, 



2t + 1 



8 m + 4nA5 k y ,k = k 



(42) 



The only scale dependence thus comes from the surface of the k- 
shell, and thus drops as k~ 2 , as explained in section |4~2l 

In the k t k' case, we find that in the non-linear regime, the 
I = point is at least two orders of magnitude above the other even 
£, and 18 orders above the odd £. The results are presented in the 
middle part of Fig. [9] for k^j ~ 1.0/iMpc -1 . The error bars were 
obtained from a bootstrap resampling of the angular dependence of 
the covariance matrix, measured from 200 Poisson densities. When 
k = k', we find that the zero-lag point effectively shifts the whole 
distribution upwards by an amount equivalent to AnC Po " ssm (k, k, 0). 

Finally, we compare the Cc distribution measured from Gaus- 
sian fields to an analytical prediction, obtained from rEqs. 1411341 . 



Cf^'(.k,k') = 2n^^(l + (-iy)S w 



(43) 



which is null for odd-^". 

We measured the Gaussian Cc from the covariance matrix of 
200 Gaussian random fields, as outlined in section |4~3l We show 
the results in the bottom part of Fig. [9] for the case where there is 



6.3 Measuring C[(kj, kj) from simulations 

We now present in this section the multipole decomposition of the 
C(k h kj, 8) matrix measured from our simulations. 

We show in Fig.[l0]the first few non-vanishing multipole mo- 
ments (i.e. t = 0, 2, 4, 6), in the case where both scales are ex- 
actly equal. We observe that higher multipoles become closer to 
the Gaussian prediction given by [Eq. 1431 , and in fact only the 
first three differ enough to have a non-negligible impact. All the er- 
ror bars in this section were obtained from bootstrap re-sampling. 
As we progress deeper in the non-linear regime, we expect to en- 
counter a mixture of the following two effects: an increase in the 
number of t required to specify the Q distribution, or in the de- 
parture from the Gaussian predictions of the Cc at fixed I. As seen 
from Fig. [Toj the departure between the multipoles and the Gaus- 
sian power increases for higher £-modes, an effect prominent in 
the first multipole. The departure becomes more modest for higher 
multipoles, and eventually we cannot distinguish between Gaussian 
and non-Gaussian. This suggests that the non-Gaussianities are en- 
capsulated in the second of the effect above mentioned. 

We then show in Fig. Qj] the same multipole moments, this 
time for the case where one scale is fixed at k = 0.61/jMpc -1 , while 
the other is allowed to vary. Once again, higher multipoles have 
smaller amplitudes, and approach the Gaussian prediction off zero. 

On the diagonal, the relative difference between the multipoles 
in the linear regime becomes smaller and converge to the predicted 
value, as expected. In addition, in the linear regime, the angular 
power of the off-diagonal elements (i.e. fc, t kj) is one to two or- 
ders of magnitude smaller than the diagonal counter part. As we 
progress to the non-linear regime however, the off-diagonal ele- 
ments decrease less rapidly, and a convenient way to express this is 
to look at the cross-correlation coefficient matrices. 



6.4 C c (k,k') matrices 

As mentioned above, we need to look at an individual I, for all 
scales combinations (k,k r ), and find the multipole beyond which 
the off-diagonal elements become negligible. The whole purpose 
behind this is to model the full covariance matrix as: 



C(k,k',6) = — Y(2C+ l)Cc(k,k')P c (cos8) 

(=0 



(44) 



where the lower I terms are measured from our simulations, and the 
others obtained from the Gaussian analytical prediction ([Eq|43)). 

In the figures of this section, we present these 'C/ matrices, 
normalized to unity on the diagonal. These are thus in some sense 
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is the Gaussian prediction, obtained from lEq. |43| . From this figure, we 
observe that the diagonal of multipoles higher than t = 4 converge to the 
Gaussian predictions. 
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Figure 12. Co matrix, normalized such that the diagonal elements are 
equal to unity. This matrix is completely equivalent cross-correlation co- 
efficient matrix of angle averaged P(k), up to a factor of 4n. It represents 
the correlation between different scales, and shows that scales smaller than 
k ~ 1.0/iMpc are correlated by more than 80 per cent. 
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Figure 11. Same as Fig.[l0] but with kj = 0.61/iMpc~' being held. The 
Gaussian prediction is zero in this case. The measurements are normalized 
by the square root of their diagonal contributions, such as to show the rela- 
tive importance of each multipole. As I increases, the off-diagonal contribu- 
tion becomes smaller, even for combinations of scales similar in amplitudes. 
The fourth point starting from the left is identical to unity for all multipoles, 
as it corresponds to a diagonal point. 



equivalent to cross-correlation coefficient matrices. Fig.ll2|presents 
the normalized Co matrix, which shows a structure similar to that of 
the cross-correlation coefficient matrices obtained previously (Fig. 
[TJ. The resemblance is not surprising, since Co = 4nC(k,k'), (we 
can see this by setting P [=0 = 1 in fEq. 1401 ) The Co(k,k') matrix 
thus contains the information about the error bars of angle averaged 
power spectra, as well as their correlation. 

By looking at the fractional error between the Co matrix and 
the actual covariance matrix of angle averaged power spectra (Fig. 
113b . we find that our method provides a very good agreement in the 
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Figure 13. Fractional error between the Co matrix obtained from the 
C(k,k',6) for all shell combinations, and that obtained directly from the 
angle averaged P(k). We do not show the largest scales, which are noisy 
due to low statistics. We have also divided the Co matrix by (4rr) for the two 
objects to match exactly. 



trans- and non-linear regimes, down to the few percent level. We do 
not show the largest scales, in which our method is more noisy, for 
reasons already explained. We recall that an extra contribution to 
Co(k, k'), not included here, comes from the non-Gaussian Poisson 
uncertainty, as discussed in section |4!2l and needs to be included in 
the final analysis. 

We now present the next few multipole matrices, and we found 
that beyond t = 4, very little information contained in the off- 
diagonal elements. Fig. Q4] shows the C2 matrix, again normalized 
to the diagonal for visual purposes. We observe that the smallest 
scales are correlated up to 60 per cent. As discussed in sections [84l 
and[9] this matrix is a requirement for an accurate treatment of BAO 
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Figure 14. C2 matrix, normalized such that the diagonal elements are equal 
to unity. The off-diagonal elements are still correlated at least at 40% for 
scales smaller than k = l.OftMpc -1 . 
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Figure 15. C4 matrix, normalized such that the diagonal elements are equal 
to unity. The off-diagonal elements are correlated at the 30% level in the 
non-linear regime, and not far off the diagonal. 



and weak lensing data analyses, first because the uncertainties on 
the underlying galaxy surveys are not isothropic, and second be- 
cause th e optim al window function is often also angle-dependent 
lLuetal.l ( f2oToh . 

Fig.ll5lshows that the correlation in the C4 matrix is still of the 
order 50 per cent for a good portion of the non-linear regime. The 
new feature here is that the strenght of the correlation of strongly 
non-linear modes among themselves starts to decrease as we move 
away from the diagonal. Fig. [T6] shows that in Ce, the matrix is 
mostly diagonal. As we progress through higher multipole mo- 
ments, the off-diagonals become even dimmer, hence do not con- 
tain significant amount of new information. From this perspective, 
a multipole expansion up to £ = 4 is probably as far as one needs 
to push in order to medel correctly the off-diagonal elements. 

Following [Eq(44), we thus propose to reconstruct the full 
C(k,k',9) from a combination of a) fully non-linear C((k,k') ma- 
trices (for { < 4), presented above, b) analytical terms given by 
[Eq. 1431 (which we scale up by 30 per cent as mentioned in sec- 
tion ^, 3t . and c) non-Gaussian Poisson error, which depends solely 
on the number density od the sampled fields. In the next section, 
we decompose and simplify these Q matrices into a handfull of 
fitting functions, and show how one can easilly reconstruct the full 
C(k, k', S) at the percent level precision. 

We next present in Fig. [T7] the ratio of the diagonal of these 
matrices to the Gaussian prediction. We observe that all of them are 
consistent with the Gaussian prediction in the linear regime. As we 
progress towards the non-linear regime, the largest departure comes 
from the Co matrix, by a factor of about 40 near k = l.O^Mpc -1 . 
We observe a turn over at smaller scales, which is caused by our 
resolution limit. We opted not to model it in our fitting formula. C2 
and C4 mildly break away from Gaussianity by factors of 4 and 2 
at the same scale. All the higher Cs are consistent with Gaussian 
statistics. Over-plotted on the figure are fitting formulas, which are 
summarized in TableQ] 
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Figure 16. matrix, normalized such that the diagonal elements are equal 
to unity. We observe that the matrix is mostly diagonal, and thus decide to 
treat C$ and all higher multipoles as purely Gaussian. 



Table 1. Fitting formulas for the ratio between the diagonals of the C((k, k) 
and the Gaussian prediction. For all £'s, the function is modeled by V(x) = 
l.0 + (x/af. 
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0.2095 


1.9980 
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0.5481 


0.7224 
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1.6025 


1.0674 
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Figure 17. Ratio of the diagonal elements of a few Cf matrices, compared 
to the Gaussian prediction. The error bars were obtained from bootstrap 
re-sampling. Over-plotted are fitting functions, summarized in Table[T] 
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Figure 18. Fractional error between the original Co matrix and that pro- 
duced with the principal Eigenvector. We do not plot the largest scales, 
which are noisy due to low statistics. 



7 FACTORIZATION OF THE C c MATRICES 



In this section, we simplify even further our results with an Eigen- 
value decomposition of the C{(k,k') matrices, normalized to unity 
on their diagonal, as shown in the figures of section l6!4l We perform 
an iterative process to factorize each matrix into a purely diagonal 
component and a symmetric, relatively smooth off-diagonal part. 
The later can be further decomposed into a small set of Eigenvec- 
tors U A (k), corresponding to the largest Eigenvalues A. These are 
then fitted with simple formulas. Combined with Gaussian predic- 
tions and fitting formulas from the deviation on the diagonal, one 
can fully reconstruct each of the C[{k, k') matrix, and thus recover 
C(A,Jf,0)as well. 

We start off the iteration by assigning the identity matrix to the 
diagonal component, which we subtract from the original matrix. 
We then extract from the remainder the principal Eigenvectors and 
recompose an new matrix as 

C ( (k,k') 



r ( {k,k') = 



iC[(k,k)C((k',k>) 



— = 8 W + V AU A {k)U A {k') 



(45) 



For the next iteration, we model the diagonal as S^—^x A(U A (k)) 2 , 
and decompose the remainder once again. We iterate until the re- 
sults converge, which takes about 4 steps. We vary the number of 
Eigenvalues we keep in our iteration, and keep the minimal num- 
ber for which the reconstruction converges. In the end, the rc(k, k') 
matrix is modeled as: 



r t {k, if) = 8 kk> [\ - AUj(k)] + ^ AU A (k)U A (k') 



(46) 



We show in Fig.[T8]the fractional error between the original matrix 
and the factorized one. The factorization of the Co matrix with one 
Eigenvector reproduces the original matrix at the few percent level. 
The same procedure is also applied for the higher multipoles, in 
which we have included the first four Eigenmodes, and we find 
that the fractional error between the reconstructed and the original 
matrix are also of the order of a few percent. 

We next fit these Eigenvectors with simple functions: for all 
£ 's, the first Eigenvector is parametrized as (J (k) = akP(y - 6k), and 
all the other vectors as U(k) = a^^mic/k 5 ). The values of (a,yS, y, 6) 
for the lowest three f s are presented in Table[2] We required that all 
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Figure 19. Fractional error between the original Co matrix and that pro- 
duced with the fitting formulas. We do not show the largest scales, which 
are noisy due to low statistics. 



of these formulas vanish as k — » 0, since the Q matrices become 
diagonal in the linear regime. The Eigenvectors of the C4 matrix 
are presented in Fig. [20] over-plotted are the fitting formulas. The 
pixel-by-pixel agreement between the original matrices and those 
obtained from the fitted formulas is with less than 10 percent for all 
t < 4 and k > 0.5. 

Larger scales fluctuate much more as they are less accu- 
rately measured, hence the pixel-by-pixel agreement is not ex- 
pected there. In addition, the matrices with £ > 6 are much harder 
to express with a small set of Eigenvectors, since the Eigenvalues 
are not decreasing fast enough. In any case, the first three harmon- 
ics we provide here contain most likely all the information one will 
ever use in realistic surveys and forecasts. 
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Figure 20. Four principal Eigenvectors of the C4 normalized matrix (solid 
lines), and corresponding fitting formulas (dotted lines). 



Table 3. Fitting parameters for the diagonal of the C p " issm (k, k') matrix, and 
for the fit of the Eigenvectors of the corresponding cross-correlation coeffi- 
cient matrix. For all three number densities (i.e. «i,2,3 = 5.0 X 10~ 5 , 1.52 X 
1(T 4 and 1.0 X 10~ 2 respectively), the Eigenvector is parametrized as 

U P f°™ son (k) = a(j +yj , and the ratio of the diagonal to the Gaussian 
prediction is fitted with y Pl >' ss '>"(k) = e e &°\ Top to bottom corresponds to 
increasing density. 
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jS 
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S 


e 


cr 


52.02 


1.0193 


0.0947 


2.1021 


2.5861 


2.6936 


2.1347 


45.09 


0.9987 


0.2034 


2.1553 


2.3407 


1.6533 


2.1965 


24.41 


0.2966 


3.3736 


0.6099 


0.6255 


-0.4321 


2.0347 



P (k) 

where C^ s °"(k,k) = ^kml \ The best-fitting parameters are 
summarized in Table [3] and the performance of the Eigenvector 
fit can be found in the Appendix. 



Table 2. Fitting parameters for the Eigenvectors of the Cc, with their cor- 
responding Eigenvalues. For all Cs, the first Eigenvector is parametrized as 
U(k) = a ( j + yj , and all the other vectors as U (k) = akPsm(yk s ). These 
were obtained from dark matter N-body simulations, but the method is gen- 
eral, and a different prescription of galaxy population may result in slightly 
different values. 
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61.9058 


0.0501 


0.0207 


0.6614 


2.3045 


2 


35.7400 


0.273 


0.8266 


1.962 


0.816 




4.4144 


0.15772 


2.4207 


0.79153 


0.032207 




1.7198 


0.14426 


4.0613 


0.76611 


-0.26272 




0.9997 


0.14414 


5.422 


0.84826 


0.31324 


4 


22.0881 


0.060399 


0.10344 


0.64008 


2.2584 




4.5984 


0.1553 


2.3370 


0.9307 


-0.1154 




2.2025 


0.1569 


3.6937 


0.92304 


0.04006 




1.4062 


0.15233 


5.1617 


0.8899 


-0.14503 



7.1 Non-Gaussian Poisson noise 

The non-Gaussian Poisson uncertainty, whose construction was 
presented in section l4~2l can conveniently be incorporated in an 
analysis by finding the corresponding Eigenvalue and Eigenvectors 
of C Po ' sso "(k,k'). Higher multipoles are not relevant as the angular 
distribution is flat, as shown in the middle plot of Fig. [2] We tested 
three number densities, corresponding to n = 5.0x 1(T 5 , 1.52x 1CT 4 
and 1.0 X l(T 2 /i 3 /Mpc 3 . In all cases, we computed the covari- 
ance matrix, decomposed into a diagonal component and a cross- 
correlation coefficient matrix, found the matrix's leading eigen- 
value and Eigenvector, then fitted the latter with the same fitting 

function: U p ° t ' sso "(k) = a(j +yj . The diagonal was also fitted 
with a simple power law of the form 



y Poisson Q\ . 



"(k,k) 



7.2 Recipe 

Here we summarize our method to generate accurate non-Gaussian 
matrices. The full C{k, k',6) matrix is then given by lEq. 1441 . where 
the I < 4 terms are obtained from the fitting functions, and the 
higher multipole moments are obtained directly from rEq. 1431 . The 
sum over these Gaussian terms can be evaluated analytically as 

i j^ + l)(l + (-l)VHj") = tf fl (i+/i) + (&,(i-/i)- 

1 - 5P 2 (m) - 9P 4 ( M ) (48) 

For the non-Gaussian terms, we proceed as follow: Each of the 
Cc(k,k') matrices, normalized on the diagonal, can be constructed 
from the first set of fit functions U,\{k) provided in Table [2] by fol- 
lowing lEq. 1461 . The 'un-normalized' C{(k,k') terms are then con- 
structed by inverting rEq. 1321 , where the diagonal elements are ob- 
tained from the product of the V((k), also summarized in Table [2] 
with the Gaussian prediction rEq. 1431 . In other words: 

C e (k,k') = (s kk {l-J^AUl e (k)y^AU A Ak)U A Ak')y 

yjv e (k)Vf(k')C° au!!S (k)Cf m,ss (k') (49) 
The complete covariance matrix is given by: 

C(k,k',fi) = — Ym+l)C { (k,k!)PM+ 

m^-(s D (l + S D (1 -ix) - 1 - 5P 2 (p.) - 9/^)) (50) 

with yu = cos(#). This can be written in a more compact form as 
C(k,k?,n) = C Gnl ,„(fc)<5(k-k') + 

^{2e+\){G t (k)S(k-k') + H ( {k,k')P f ( i i^ (51) 



^Poisson 
Gauss 



(k,k) 



e € k" 



(47) 



with 

G c (k) = C GauS s(k){V e {k) - 1) 



(52) 
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and 



n 



F AX (k)F A Ak')-Fl e (k)6(k-k 



'>) 



FaAV = U A Ak)^AV ( (k)C Gm , s (k) 



(53) 



(54) 



We conclude this section with a word of caution when us- 
ing the fitting formulas provided here, in the sense that the range 
of validity of the fit has not been tested on other cosmological 
volumes. Consequently, we advice that one should limits itself to 
k < 2.0AMpc _I . 



8 MEASURING THE IMPACT WITH SELECTION 
FUNCTIONS 

This section serves as a toy model for a realistic non-Gaussian er- 
ror analysis, as it incorporates the non-Gaussian covariance ma- 
trix measure d from N-body simu lations with the 2dFGRS selec- 
tion function (Norberg et al. 2002). We compare the estimated error 
bars on P(k) between the naive, purely diagonal, Gaussian covari- 
ance matrix, the effect of the one dimensional window function as 
prescribed by the FKP formalism, the unconvolved non-Gaussian 
covariance as measured from our 200 N-body simulations, and the 
convolved non-Gaussian matrix. 

We recall that in a periodic volume, a selection function that 
is exactly uniform throughout the volumes makes the observed and 
true covariance matrices exactly equal. That is only true in simu- 
lated volumes, and in that case, the optimal estimator for the power 
spectrum covariance matrix is obtained from the non-Gaussian co- 
variance matrix directly measured from the simulations. 

Non-periodicity is best dealt with by zero-padding the ob- 
served survey, which results in some coupling between different 
measure power spectrum bands. The coupling becomes more im- 
portant as the selection function departs from a top hat, and in that 
case, the best estimator of the observed covariance matrix is a con- 
volution of the 6-dimensional covariance over both vectors (k, k'), 
given by: 

_ _ £ k » u» C, TO (k", k"')|W(k - k")| 2 |tV(k' - k"')| 2 

Cofcs(k ' k } " (^7v c ZxW 2 (x)wW { ' 

The denominator is straight forward to calculate, while the numer- 
ator is a 6-dimensional integral, which must be calculated at all of 
the 6-dimensional coordinates, a task is computationally impossi- 
ble to perform. For example, assuming that the size of the grid is 
n 3 cells, then for each (k, k') pair, we have to sum over n 6 terms. 
We also have n 6 such pairs, and each term takes about 3 flop. For 
n = 100, this process would take 3 * 10 24 flop, and current super- 
computers are in the regime of resolving 10 12 flop per seconds. 
The above calculation would therefore take about 3000 years to 
complete. With the factorization proposed in this work however, it 
is now possible to break down the computation into smaller pieces 
and reduce the loops to 7 dimensions at most. 



8.1 Factorization of the 6-dimensional covariance matrix 

We propose here a break down of the true covariance matrix 
C(k",k"') into a product of simple functions of the form H t (k"), 
Gc(k") and P{(ji) where the angular components come exclusively 
from the Legendre polynomials. As explained in section [2~3l the 
argument of the Legendre polynomials, fi, is the (cosine of ) an- 
gle between k" and k'", which must first be expressed in terms of 



Table 4. List of weights w(9, <j>) needed for the angular integrals over the 
selection function. These can be precomputed to speed up the convolution. 
All integrals are in the form of rEq |57| . 



cos 2 (<9) sin 2 ((?) 


cos 2 (e)e ±2,t * 


sin 2 (2e)e ±! '* 


cos 4 (<9) sin 4 ((?) 


sin 4 (0)e ±2l '<* 


sin 4 (260e ±4i * 


sin 2 (20) sin 2 (20)e ±2 '* 


sin(0)cos 3 (e)e ±! '« i 


cos(e)sin 3 (e)e ±/ * 



(0",(f>",0"',(f>"'), following rEq. [T4lH The only multipoles that ap- 
pear in our equations are t = 0, 2, 4, so fi is to be expanded at most 
up to the fourth power. For a full factorization, the terms including 
cos(0" - (/>"') must further be re-casted in their exponential form 
with Euler's formula. 

When computing the convolution, the first term on the right 
hand side of fEa lS II is spherically symmetric, hence it must be 
convolved with the selection function as: 



ctLAKk') = V c Gauss (k"W(k" - k)\ 2 \w(k" - k'f 



k" 



(56) 



which is pretty much the FKP prescription, namely that the selec- 
tion function is the only cause of mode coupling. 

When computing the other terms of fEq |5 1 1 . which are the 
non-Gaussian contributions, we use the fact that the only coupling 
between the k" and k'" vectors comes from the delta function, 
which couples their radial components. This means that all the an- 
gular integrations can be precomputed and stored in memory. For 
example, the only angular dependence in the C = multipole comes 
from the selection function itself, hence we can precompute 



X(k,k") = |W(k-k")| 2 sin(0")w(6»",0") 



(57) 



6" 4" 



and the convolution is now four dimension smaller. The weight 
function w(d", 0") is equal to unity for the Co term, and the sin(#") 
comes in from the Jacobian in angular integration. For the higher 
multipoles, more of such terms must be precomputed as well, 
whose weight functions are summarized in Table|4] 



8.2 The 2dFGRS selection function 

The 2dFGRS dColless et al.ll2003h is comprised of two major re- 
gions, the NGP and the SGP, each of which take the overall form 
of a fan of 75 x 5 degrees, extending from z = 0.02 to 0.22. The 
selection function is constructed by first integrating the luminosity 
function d<b(L)/dL over all the observed luminosity range, which is 
both redshift and angle dependent, and multiplying the result by the 
redshift completion function R(9, (f>). Namely, we define the galaxy 
number density h(z, 0, cf>) as: 



n(z,6,tp) 



"\i , „„-„(;,fl,i« 



dL 



ydy 



(58) 



where y = LJL*, L* being the characteristic galaxy luminosity, and 
-gf- given by 



d<SKy) 
dL 



(59) 



5 In this section, we use fi instead of cos($) to denoted the (cosine of the) 
angle between the two Fourier modes, to avoid confusion with ff and ff" , 
which corresponds to the angle of (k", k'") with respect to the jc-axis. 
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The three parameters C>*, a and M* - 5log 10 h are obtained from 
the 2dFGRS as -1.21, 1.61x10^ /i 3 Mpc~ 3 and -19.66 respectively. 



The integral over y gives an incomplete gamma function: 

n(z, 9, 4>) = <S>*T(a + 2, y, nm (z, 9, </>)) 

in which the term y mj „ can be expressed as 

lOgioOminfo e > 0)) 



(60) 



0.4(m* -51og 10 /t- 6,(0,0)+ 
/ D L {z) \ z + 6z 2 \ 

51ogl Hl^J + TT20?) 



(61) 



where b,{9, (f>) is the angular dependence of the magnitude sensitiv- 
ity, D L (z) is the luminosity distance that is used to convert between 
absolute relative magnitudes, and the last term in the right hand 
side is the X"-and e-corrections. Finally, the selection function is 
simply : 



W(z, 9,<p) = h(z,9,(f>)R(9,(l>) 



(62) 



Both R{9,<f>) and bj(9,(f>) are publicly available from the 2dFGRS 
websit^]. It is possible to obtain an even more accurate selection 
function by taking into account the redshift dependence of the mag- 
nitude sensitivity, however we do not need such an accuracy for the 
current work. We finally normalize the selection function such that 



\W(k)\ 2 d 3 k = 1 



(63) 



To understand the impact of the non-Gaussian Poisson un- 
certainty on the measured uncertainty, we test various templates, 
keeping the 2dFGRS selection function fixed. We follow the pro- 
cedure of section l4~2l with an average number density of n ga i = 
1.52 x l(T 4 /! 3 Mpc~ 3 , which corresponds to an early data release 
of the 2dFGRS data. The final release contains more objects, and 
has a density of about n = 5.0 x 10~ 2 /i 3 Mpc~ 3 . By comparison, 
the Poisson uncertainty corresponding to the number count of the 
Wiggle-Z survey could be modeled with n = 5.0 x 10~ 5 /i 3 Mpc~ 3 for 
partial data and about 2.0 x 10~ 4 /? 3 Mpc~ 3 for the final data release. 
We thus opt for two more number densities: n = 1.52 x 10~ 4 and 
n = 1.0 x 10 2 . 



8.3 Results 

We assign the selection function on to a 256x256x128 grid, where 
the lower resolution is along the direction perpendicular from the 
NGP. We precompute the Fourier transform, W(k) and square each 
terms. Fig. [21] shows a comparison between the angle average of 
|W(k)| 2 and a fitting function provided by the 2dFGRS group. 

We then define a second set of bins in spherical coordinates, 
over which we perform the convolution. For that purpose, we divide 
the original volume of the survey into 64 radial bins, 48 polar bins 
and 32 azimuthal bins. The selection function is assigned in the 
grid by averaging over the 27 closest cells in the original grid. We 
have included a sin(S) terms in each integrals over the polar angle, 
and a k 2 in each radial integral to properly account for the Jacobian 
matrix in spherical coordinates. 

Fig. [22] shows the diagonal of the convolved covariance ma- 
trix, divided by P 2 (k), for the FKP prescription and for the progres- 
sive inclusion of £ = 0, 2 and 4 multipoles. Also overploted is the 
non-Gaussian results without the convolution. We see that already 
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Figure 21. The angle average of the 2dFGRS selection function, compared 
to an approximate fit provided by jPercival et alj200ll) . The fit is not perfect 
as it was obtained with an earlier estimate of the selection function. We also 
note that our method differs in details with that used in <Cole et alj : 2005) 
by that fact that we imposed a cut at redshift of z = 0.22, and that we used 
a somewhat smaller resolution. 



at k ~ 0.1/zMpc~', the non-Gaussian fractional error, after the con- 
volution, deviates from the FKP prescription by a factor of about 
3.0, while the unconvolved Co still traces quite well the FKP curve. 
This means that the mode mixing caused by the convolution with 
the survey selection function increases significantly the variance of 
the observed power spectrum. The departure gets amplified as one 
progresses towards higher k— modes, and by k ~ 1.0, the uncon- 
volved Co departs from the FKP prescription by almost two orders 
of magnitudes. Interestingly, the convolved Co merges with the un- 
convolved counterpart at k ~ 0.5, where the BAO scale is usually 
cut off. inclusion of higher multipole increases the variance by a 
factor of about 2.0. 

We have overplotted a simple smooth fitting function of the 
form : 

2.3 

Cf it {k) = CAk)(l + r-, r-r + 0.0007) (64) 

1 y gy (0.08/yt) 3 - 7 + (0.08/fc) 1 1 

which approximates the contribution from the three lower multi- 
poles. 

Fig.[24]shows the convolved cross-correlation coefficient ma- 
trix, where the angle average has been taken after the convolution. 
It is also possible to factorize this matrix, hence we proceed to an 
Eigenvalue decomposition, following the same iterative procedure 
as in section [7] solving for the first Eigenvector only. The Eigen- 
value was found to be A = 19.7833, and we used the sum of a 
quadratic and a Gaussian function to model the Eigenvector: 



UT(k) 



Aexp[ — jlog 2 (k/k„)]+ 
(alog 2 (k/k ) + blog(k/k ) + c) 



(65) 



6 www.mso.anu.edu.au/2dFGRS/ 



with A = 0.1233,o- = 1.299, a = 0.0049, b = 0.0042, c = 0.0052 
and (k p ,k ) = (0.17, 0.008)/!Mpc~' respectively. A comparison of 
the fit and the actual vector is presented in Fig. [23] The noise re- 
duced cross-correlation coefficient matrix is presented in Fig. [25] 
We observe that the Fourier modes are already more than 50 per 
cent correlated at k = 0.1/iMpc~', an significant enhancement com- 
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Figure 22. Diagonal of the convoluted covariance matrix, first with no 
multipole i.e. following FKP prescription (bottom-most dotted line), then 
with the progressive inclusion of the Co (solid black points), the C2 (open 
circles) and the C4 multipoles (thick solid line). Also shown is the diagonal 
of the unconvolved Co terms directly measured from N-body simulations 
(dashed line), and a fitting function for the total covariance (thin solid line). 
Finally, the inclusion of the non-Gaussian Poisson noise is represented by 
three dotted lines, representing the three number density detailed in Table 
[5] The 2dFGRS final data release has a number density of the order 5.0 X 
10~ 2 /i 3 Mpc~ 3 , which thus lies between TI2 and 113. 
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Figure 23. Principal Eigenvector of the convolved Co matrix, compared 
to a simple fitting formula. The fractional error of the fitting function is at 
most 13 per cent. 



pared to the unconvolved Co matrix, in which the equivalent cou- 
pling occurs roughly towards k = 0.22/iMpc -1 . This would most 
likely have an impact on a non-Gaussian BAO analysis. 



8.4 Applications to weak tensing 



The results presented in section 16.41 and the recipe presented in 
the previous section can find useful applications in the field of 
weak lensing. Convergence maps, for instance, are constructed 
from a redshift integral over a past line cone filled with dark matter, 
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Figure 24. Normalized convoluted covariance matrix with all three multi- 
pole. 
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Figure 25. Normalized convoluted covariance matrix with all three multi- 
pole, reconstructed from a fit of the principal Eigenvector. 



weighted by a geometric kernel. Because of the projection nature of 
this process, the survey maps are sensitive to both large and small 
scales, where non-Gau ssianities have be en observed in the conver- 
gence power spectrum foore et al.ll2 009). 

It has recently been demonstrated that weak lensing of high 
redshift sources by large scale st ructures can se rve as a probe of 
the dark energy equation of state (Huterei 2002). To make a com- 
plex story short, the weak lensing power spectrum is closely related 
to that of the three dim ensional dark matter via Limber's approxi- 
mation jLimbed[l954l) , hence it is similarly sensitive to cosmolog- 
ical parameters. It was then realized that the same structures were 
also distorting the signal of high redshift source of 21 cm emis- 
sion dMetcal f & White 2009), hence the detection of weak lensing 
from diffuse temperature fields could also constrain the dark en- 
ergy. The lensing fields are quadratic functions of smoothed tem- 
perature fields, and the optimal smoothing window function de- 
pends not only on the the parameter und er study, but als o on the 
statistical nature of the source and lenses Jl"u & Pen 2008). 
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Optimal quadratic estimators of len sing fields were first 
obtained under the Gaus sian assumption (Hu & Oka motol 1 20021 ; 
IZahn & Zaldarriagai r2006). but it was soon realized that when the 
observed field is non-Gaussian, the proposed estimators are no 
longer optimal: the noise estimation can be underestimated by sev- 
eral orders of magnitude l lLu & Penll2008l) . Optimal non-Gaussian 
lensing estimat ors were then obtained from N-body simulations 
jLu et alj|2010h , and it was found that the optimal smoothing win- 
dow function for dark energy involves the first two multipoles of 
the dark matter power spectru m covariance m atrix, Co(k,k') and 
C 2 (k,k') (see [Eqs. 23 - 24] in dLu et alj|201oh ). The tools devel- 
oped in the present paper thus allow one to construct, for the first 
time and from simple fitting functions, optimal non-Gaussian esti- 
mators of dark energy parameters from 21 cm temperature maps. 

The survey selection function was factored out of the above 
discussion, which was interested solely in the non-Gaussianities in- 
herent in the sources and lenses. However, cosmology from weak 
lensing observations of galaxy and quasar surveys is affected by 
selection functions, in a similar way to BAO analyses discussed in 
this paper: the underlying 2- and 4-point functions in Fourier space 
are inevitably convolved with the survey selection function. With 
the factorization proposed in the paper, this convolution should also 
be relatively simple to perform, and allows one to measure non- 
Gaussian error bars on the lensing power spectrum, which includes 
the effect of the non-linear dynamics, of the selection function, and 
possibly of Poisson noise. 

In particular, the weak lensing convergence field k is obtained 
from the dark matter field 5 from 

« = f Mx>X.tWx (66) 
Jo 

with 

3H?:Q m I y \ 
Mx,Xs) = -£s~X 1 1 - "J d + zfc» ( 67 > 

H is the Hubble parameter, c the speed of light, x and Xs the co- 
moving distances to the lenses and to the source respectively. When 
one has knowledge of the three dimensional selection function, it 
can be incorporated in the above equation as 

k«b,4>)= ['wQc,x.Wx, 0,4>)W(x,e,(/>)d x (68) 

Jo 

we can be further compactified if one absorbs w in the definition 
of the selection function. The result takes the form of an inte- 
gral over one dimension of an observed density fields: K(d,ip) = 
f 6(x)W(x)dx. This means that all the non-Gaussian calculations 
performed on the three dimensional density fields (i.e. see section 
12.3b can be extended to weak lensing maps by the simple modifica- 
tion of the selection function mentioned above, plus an integration 
over the third dimension (in the small angle approximation at least). 
In particular, the best estimate of the non-Gaussian lensing covari- 
ance matrix is obtained from solving [Eq. 1131 with the modified 
selection function, then by integrating the two k : -components with, 
at each integration step, the inclusion of an extra weight X f° r the 
conversion of A>modes to f-modes. 



9 DISCUSSION 

We have found that even for modes of k ~ 0. lAMpc -1 , the non- 
Gaussian error bars are higher than those prescribed by the FKP 
method by a factor of a few, due to mode coupling caused by the 
convolution of the selection function. This has to be put in contrast 



with results from pure N-body simulations, which show that the de- 
parture from Gaussianity reaches this sort of amplitudes at higher 
£-modes, as seen from Fig. [22] We also observe that with the 2dF- 
GRS, the non-Gaussian Poisson noise plays an important role if the 
number density is smaller than 0.01/i 3 Mpc~ 3 , but is not enough to 
characterize all of the non-Gaussian features of the density field. 
The Co term is the leading contribution of the enhancement ob- 
served in the range k = 0.06 — 0.4/?Mpc~' , but for larger fc-modes, 
C2 and C4 both play an important role. 

Without the convolution, keeping only the Co term, and 
assuming that the BAO measurement was performed with a 
non-Gaussian estimator, the propagation of the non-Gaussianities 
on to the BAO dilatio n scale produces very similar constraints 
( Taka hashietafll201ll) . The estimators that are used in the data 
analyses however are Gaussian, while the power spectrum covari- 
ance matrices that enter the calculations are either Gau ssian or ob- 
tained with mock catalogs. As pointed out previously (Ng an et all 
1201 lh . the estimators constructed in such a way are inconsistent and 
should be noise weighted. When correcting for that effect, it was 
found that the consistent - but suboptimal - error bars are about 10 
per cent higher than those obtained assuming an optimal estimator. 
In the light of the current results, the observed (i.e. convolved) co- 
variance matrix is even less Gaussian, and it is not obvious that the 
error on the BAO dilation scale will be unaffected by this new esti- 
mates, since our measurement show significant deviations at scales 
as large as 0.1/iMpc -1 . 

It is worth mentioning again that the measurement of the 
Co(k, k') from C(k,k',8) provides an alternative way to extract the 
covariance matrix of the angle average power spectra. Although 
the mean value of both methods is identical, i.e. unbiased, the sec- 
ond gives us a better handle on the error on each matrix element, 
hence provide an optimal measurement of their uncertainty. We 
have shown in this paper that each matrix element receives its dom- 
inant contribution from small angles, while larger angles are more 
noisy. It is thus in possible to re-weight the sum by taking this new 
information into account, and obtain more accurate error bars on 
each matrix element, compared to the current bootstrap method 
(HDP2). As mentioned in the introduction, our next objective is 
to achieve a similar accuracy with a much lower number of sim- 
ulations. This would revolutionize the field of observational cos- 
mology as the covariance matrix would be measured internally, i.e. 
directly from the data. 

The techniques presented in this paper call for extensions, as 
we did not include redshift distortions nor shot noise in our anal- 
ysis. The latter will become importan t when repeating thi s pro- 
cedure on haloes, and it was shown dNevrinck et alj 120061) that 
the Fisher information in haloes is also departing from Gaussian- 
ity. It is straight forward to perform a similar analyses with a 
quadratic halo model, where the halo density is parametrized by 
Shah(x) = Ad(x) + B6 2 (x). This involves an extra cross correlation 
between the linear and quadratic term, and leaves some room for 
the choice of A and B, and ultimately, one should work straight 
from a halo catalog. The optimal estimator should also be based on 
a cosmology independent model of the covariance matrix, hence 
one should compute how the fitting functions scale with Q.„, and u. 

As mentioned earlier, the effect of the selection functions 
is enhanced for survey geometries that are different from top- 
hats, and it would be interesting to repeat some of the BAO data 
analyses that were performed on such surveys, like the 2dFGRS, 
Wiggle-Z. The current method also applies to surveys with ir- 
regular geometries like those obtained from the Lyman -o- forest 
( McDonal d~& Eisensteinll2007l : iMcOuinn & Whitell201 lh . and we 
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are hoping that it will considered in the elaboration of these future 
analysis pipelines. 

We leave it for future work to match our results with predic- 
tions from higher order perturbation theory. We would like to verify 
that the angular dependence we observe in the covariance matrix 
is predicted by a complete 4-points function analysis, at least in 
the trans-linear regime. In addition, the extraction of non-Gaussian 
error bars from two dimensional angular clustering could also be 
performed with techniques similar to those employed here. 



10 CONCLUSION 

Estimating accurately the non-linear covariance matrix of the mat- 
ter power spectrum is essential when constraining cosmological pa- 
rameters including, but not restricted to, the dark energy equation 
of state oj. So far, many BAO analyses from galaxy surveys were 
performed under the assumption that the underlying density field 
is Gaussian, which yields a suboptimal measurement of the mean 
power spectrum (and thus of the BAO dilation scale), and, at least 
as important, the error bars are biased. 

To estimate unbiased error bars on the dilation scale is a chal- 
lenging task but can now be d one. In the simple case of periodic 
volume, it was shown recently dNganetalfeOllh that, first, an un- 
biased error bar on a suboptimal measurement of the mean could be 
obtained from the knowledge of the underlying covariance matrix. 
Second, if one did measure optimally the mean BAO dilation scale, 
then the optimal measurement of the error requires an estimate of 
the inverse of the power spectrum covariance matrix. This is much 
more challenging due to the presence of noise, even when dealing 
with simulations embedded in periodic volumes, but improves the 
constraining performance by a significant amount. 

When estimating the power spectrum and its uncertainty from 
data, the survey selection function complicates the calculations 
since the observed quantities are actually convolved with the se- 
lection function. Since the covariance matrix is not isotropic, as 
it depends on the relative angle between two Fourier modes, the 
convolution cannot be simply factored into two radial components. 
Hence we are left with a challenging six-dimensional integral to 
perform, which so far has been an unresolved problem. 

In this paper, we present a method to perform this convolution 
for an arbitrary galaxy survey selection function, that thus allows 
one to measure unbiased error bars on the matter power spectrum. 
The estimate is still suboptimal, unless one combines our tools with 
the PKL formalism, but we have nevertheless removed the bias on 
the error bar. 

From an ensemble of 200 N-body simulations, we have mea- 
sured the angular dependence of the covariance of the matter den- 
sity power spectrum. We have found that on large scales, there is 
only a weak dependence, consistent with the Gaussian aspect of the 
fields in that regime. On smaller scales, however, we have detected 
a strong signal coming from Fourier modes separated by small an- 
gles. This comes from the fact that the complex phases of these 
modes are similar, hence they tend to couple first. We next ex- 
panded the covariance C(k, k', 6) into a multipole series, and found 
that only the first three even poles were significantly different from 
Gaussian fields. 

We further decomposed these Cf(k, k') matrices into diagonal 
terms and cross-correlation coefficient matrices, from which we ex- 
tracted the principal Eigenvectors. This allowed us to break down 
the underlying six-dimensional covariance into a set of Eigenvec- 
tors, Eigenvalues and three diagonals terms. We provided simple 



fitting formulas for each of these quantities, and thus enable one to 
construct a full six-dimensional covariance matrix with an accuracy 
at the few per cent level. 

Intrinsically, non-Gaussianities introduce A' 2 matrix elements 
to be measured in N-body simulations, as opposed to N for Gaus- 
sian fields. With the proposed method, the number of parameters to 
measure is reduced to a handful, even if the survey selection func- 
tion is non-trivial. This opens up the possibility to measure non- 
Gaussianities directly from the data, which we will investigate in 
part II of our work. 

This factorization is necessary in order to estimate unbiased 
non-Gaussian error bars on a realistic galaxy survey, that must in- 
clude the effect of the survey selection function. We found that in 
the case of the 2dFGRS selection function, the non-Gaussian frac- 
tional variance at k ~ 0.1AMpc~' is larger by a factor of three com- 
pared to the estimate from the FKP prescription, and more than 
an order of magnitude at k ~ 0.4/iMpc~'. With similar techniques, 
we were able to propagate a few templates of non-Gaussian Pois- 
son error matrices into the convolution and estimate the impact on 
the measured power spectrum. We showed that with the 2dFGRS 
selection function, the non-Gaussian Poisson noise corresponding 
to a number density significantly lower than 0.1h 3 Mpc~ 3 has a 
large effect on the fractional variance at scales relevant for BAO 
analyses and should be incorporated in an unbiased analysis. The 
cross-correlation coefficient matrix of the convolved power spec- 
trum shows that the correlation propagates to larger scales in the 
convolution process, and should have a larger impact on BAO anal- 
yses for instance. We conclude by emphasizing on the fact that con- 
straints on cosmological parameters obtained from BAO analyses 
of galaxy surveys are currently significantly biased and suboptimal, 
but that both of these effects can now be dealt with in further anal- 
yses. 
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APPENDIX A: NORMALIZATION 

We present in this appendix the interpretation of the normalization 
E^(Afc) of the angular covariance (see [Eq.l26land l27l ). The shape 
of this function is plotted in Fig. lAll for a scale that corresponds to 
k = 1.00/iMpc -1 , which lies close to the transition between trans- 
linear and non-linear regime at z = 0.5. 

In fact, the normalization is simply a counting of the differ- 
ent combinations of ky that produce a given Ak. We see from the 
figure that at small angles (i.e. small Ak), there is a lot of possible 
combinations, and that the function rapidly decreases with the an- 
gle. To capture this, imagine counting the number of times we can 
embed both ends of a given vector on a given spherical shell. In the 
thin shell approximation and for a vector with non-zero length, the 
vector can be moved around an axis parallel to the vector, that also 
passes through the centre of the shell. In terms of solid geometry, 
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Figure Al. Function that appears in the normalization of the angular co- 
variance, F(ki, kj, Ak) (equation [27), plotted here versus 8, for two shells of 
kjj = l.QAMpc . This geometrical factor represents the number of com- 
binations, for each angle separation, as obtained from a given set of shells. 
The dotted curve was obtained from the numerical integration of Eq |27l and 
fits very well the results obtained from the FFTW (solid curve). Curves like 
this one are obtained numerically for each scale combination k\ j . 
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Figure CI. Principal Eigenvector of the cross-correlation coefficient matrix 
associated with the non-Gaussian Poisson noise, compare to our best-fitting 
formula. 



The factor of 2tt comes from the integral over the tp angle, and Afi 
is half the distance to the first knot. 



the vector spans the flat side of a cylinder that intersects the shell. 
The normalization is proportional to the circumference of the cylin- 
der's basis, and as the length of the vector decreases, so does the 
height of the cylinder, hence its circumference increases. In the dis- 
crete case, the shells have a finite thickness, and are constructed out 
of a grid. The normalization thus produces integer counts, which 
discretizes the subtended angles. 



APPENDIX B: LEGENDRE-GAUSS WEIGHTED 
SUMMATION 

The conversion of the integral into a sum is performed using a 
Legendre-Gauss weighted sum jAbbottll2005h . in which £ 'coloca- 
tion' knots, which we label fi k with k = l,2,..,€, are placed at the 
zeros of the Legendre polynomial P[(^i). We choose £ = 101, and 
we exclude the end points at /j = ± 1 in order to isolate the zero-lag 
contribution. The weights w* are given by: 



a-4)(.dP^ m /dn(Mk)) 2 

This Gaussian quadrature gives an exact representation of the in- 
tegral for polynomials of degree 201 or less, and provides a pretty 
good fit to most of our C(kj,kj,9). In the linear regime, the dis- 
cretization effect becomes important, and the number of angles one 
can make between the grid cells drops down as k 2 . In the case were 
fewer points are available, we choose £ = 51,21, 11 or 5 depending 
on the number of available angular bins. Once we have specified the 
knots, then, for each scale combination, we interpolate the angular 
covariance on to these knots, and then perform the weighted sum. 
As mentioned above, we always treat the zero-lag point separately 
in order to avoid interpolating its value to the nearest neighbors. We 
thus break the summation in two pieces: 

q = 2^ PeUi^dki^n^k + 2nC(ki,kj,ti = 1)A//(1 + (-l) f ) 

(B2) 



APPENDIX C: EIGENVECTOR OF THE POISSON NOISE 

This Appendix presents the Eigenvector that best describes the non- 
Gaussian Poisson noise, as discussed in section rTTI We restrict our- 
selves with the case where the number density is the highest, even 
though similar analyses can be carried for the other values of n we 
studied in this paper. We present in Fig. IC II the Eigenvector itself, 
compared to the best-fitting formula provided. We next compare the 
covariance matrix constructed from the fitting functions with the 
original, and present the fractional error in Fig. IC2I which shows 
a few per cent a greement. When compared with the predictions 
from (Cohn 2006), we observe that the overall trends are consistent: 
first, the Gaussian contribution to the error decreases as one probes 
smaller scales. Second, densities with lower n see their Gaussian 
contribution being reduced in the trans-linear regime, where the 
non-Gaussian Poisson counting becomes more important. Third, 
densities with lower n produce larger cross-correlation coefficients 
of tran-linear scales, also in accordance with the predictions. 
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